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Abstract. In order to describe unbalanced ultracold fermionic quantum gases on optical lattices in a har- 
monic trap, we investigate an attractive {U < 0) asymmetric (t| 7^ t^) Hubbard model with a Zeeman-like 
magnetic field. In view of the model's spatial inhomogeneity, we focus in this paper on the solution at 
Hartree-Fock level. The Hartree-Fock Hamiltonian is diagonalized with particular emphasis on superfluid 
phases. For the special case of spin-independent hopping we analytically determine the number of so- 
lutions of the resulting self-consistency equations and the nature of the possible ground states at weak 
coupling. We present the phase diagram of the homogeneous system and numerical results for unbalanced 
Fermi-mixtures obtained within the local density approximation. In particular, we find a fascinating shell 
structure, involving normal and superfiuid phases. For the general case of spin-dependent hopping we cal- 
culate the density of states and the possible superfiuid phases in the ground state. In particular, we find 
a new magnetized superfiuid phase. 

PACS. 03.75.Hh Static properties of condensates; thermodynamical, statistical, and structural properties 
- 39.25.-f k Atom manipulation (scanning probe microscopy, laser cooling, etc.) (see also 82.37.Gk STM and 
AFM manipulations of a single molecule in physical chemistry and chemical physics; for atom manipulation 
in nanofabrication and processing, see 81.16.Ta) - 71.10.Fd Lattice fermion models (Hubbard model, etc.) 



1 Introduction 

Experimental background Experiments on unbalanced 
ultracold fermionic quantum gases have become of signifi- 
cant interest during the last year [1,2]. While these exper- 
iments were performed on quantum gases in the absence 
of an optical lattice, it will only be a matter of time be- 
fore a lattice structure is superimposed on the trap poten- 
tial. Since the experiments without a lattice have already 
been theoretically analyzed [3,4], in this paper we do the 
next step and investigate systems with an additional op- 
tical lattice. Such lattice systems can be described with 
the help of Hubbard- type models [5-12]. Hence, since we 
are interested in superfluid phases with a possible imbal- 
ance in the pseudospins' occupation numbers we analyze 
an appropriate Hubbard model with an attractive inter- 
action {U < 0) and a Zeeman-like magnetic field. The 
symmetric and magnetic field-free version of this Hamil- 
tonian has been intensively discussed in the seminal paper 
by P. Noizieres and S. Schmitt-Rink [13]. The pseudospin- 
states in this Hamiltonian represent different hyperfine 
states of one and the same atomic species or, alterna- 
tively, fermionic atoms with different masses. Accordingly, 
we also allow for a possible asymmetry in the model, i.e., 
for a possibly spin-dependent hopping {t-^ 7^ ^i)- To deal 
with the spatial inhomogeneity, caused by the presence of 
the harmonic trap, we focus in this paper on solutions at 



Hartree-Fock level; possible extensions and improvements 
are discussed in section 6. The Hartree-Fock approxima- 
tion limits predictions to the region of not-too-strong (at- 
tractive) coupling between the fermions. In this paper, 
therefore, we restrict consideration to the moderate weak 
coupling regime. Clearly this regime can be realized ex- 
perimentally, since interactions are tunable in ultracold 
quantum gases with the help of Feshbach resonances (see 
[6,8,14]). 



Model Hamiltonian A class of models, which is well suited 
for describing ultracold quantum gases on optical lattices, 
are models with a local Hubbard interaction. Here we de- 
scribe only the basic ingredients of the Hamiltionian. De- 
tails of the appropriate model and its properties are dis- 
cussed in section 2. 

In ultracold quantum gases on optical lattices, the Hub- 
bard interaction is generated as follows. The periodic po- 
tential induced by the coherent laser beam has the form 



y(x) 



cos 



-) 



(1) 



where Vq is proportional to the laser beam's intensity, 
a = -J is the lattice constant and A is the laser beam's 

wavelength. In the tight binding limit, Vq 3> 5^^, only 
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the local overlap matrix element U of the interaction con- 
tributes to the Hamiltonian. At sufficiently low tempera- 
tures (fceT ^ ^'gap, where iJgap is the band gap between 
the lowest two Bloch bands) an effective one-band Hub- 
bard model may then be derived from the potential (1). 
Apart from the on-site interaction term, which is here as- 
sumed to be attractive, the Hubbard Hamiltonian contains 
a nearest-neighbor hopping term. The hopping amplitudes 
will in general be pseudospin-dependent, t| ^ since 
the pseudospin index labels atoms in different hyperfine 
states or with different masses (sec [15]). Furthermore, the 
particle density in the grand canonical ensemble is tuned 
with the help of a chemical potential, and the imbalance 
in the occupation numbers of the pseudospin species is 
regulated by a Zeeman-like magnetic field. 



Structure of the paper The plan of the paper is as fol- 
lows. First, in section 2, we describe the Hamiltionian, 
its symmetry properties and the relevant self-consistency 
equations for the particle density and the superfluid or- 
der parameter. Then, in sections 3 and 4, we discuss the 
possible solutions of the self-consistency relations for sym- 
metric hopping, t-f = ti, and the corresponding numerical 
results. Our most important result here is the finding of 
an interesting shell structure for the spatial distribution of 
the quantum gas in the trap. Results for the general case 
of asymmetric hopping, t| ^ t^, are discussed in section 5. 
Here we find a new magnetized superfluid phase. We end 
(in section 6) with a brief summary and an outlook. 



2.1 Symmetry properties 

Chemical potential at half filling. The discrete symme- 
tries of the standard Hubbard model (B = and t| = t^) 
are well known from the literature. In this case, the grand 
canonical Hamiltonian K is invariant, e.g., under a spin- 
exchange transformation and. at half filling on 

O l(T 1, — (7 ' O 

bipartite lattices, also under a general particle-hole trans- 
formation c|g. — > (— 1)' Cjg., where ( — 1)' equals -1-1 for sites 
i on the A-sublattice and —1 for i on the B-sublattice. 
From these discrete symmetries it can be derived that, at 
half filling, the chemical potential is given by 



In the general case (i.e., both B ^ and t-^ 7^ ij,) these 
symmetries do not exist, so that at half filling the chemical 
potential y^t is a more complicated function of the system 
parameters. However, if either B = Q or = ^i: equation 
(3) is still vahd: For B = (but possibly ^ ti), the 
Hamiltonian is not invariant under a spin-exchange trans- 
formation, but it is invariant under a general particle-hole 
transformation at half filling. Hence, (3) is still valid. If 
t| ~ ti (but possibly B ^ Q) the Hamiltonian (2) is nei- 
ther invariant under a general particle-hole transformation 
nor under a spin-exchange transformation; nevertheless it 
is mapped onto itself at half filling by consecutively per- 
forming both transformations, so that (3) is also satisfied 
in this case. 



2 Properties of the Hamiltonian 

The full grand-canonical Hamiltonian K ~ H — fiN, de- 
scribing (possibly spin-dependent) hopping of two differ- 
ent fermionic atomic species in a (possibly unbalanced) 
mixture, reads: 

(ij)- i (2) 

- ^ ^ rii^ -I- B ^ (JTlia ■ 

1(7 1(7 

The various terms represent the hopping of the fermions, 
their local Hubbard interaction, the chemical potential 
and a Zeeman-like magnetic field. Alternatively, the last 
two terms can of course also be interpreted as a spin- 
dependent chemical potential /icr = /i — aB. The summa- 
tion label (ij) in the kinetic energy indicates that all near- 
est neighbor sites i and j are summed over, so that every 
bond (ij) is counted twice. Throughout, we assume that 
the lattice has (hyper)cubical structure, although many 
symmetry properties are more generally valid for bipartite 
lattices. In the following, we first discuss some symme- 
try properties of the Hamiltonian (2), and then we derive 
the self-consistency equations within the Hartree-Fock ap- 
proximation. 



Relation to the repulsive U model. In order to under- 
stand the properties of the model (2) at J7 < and ar- 
bitrary filling, it is often helpful to map the Hamiltonian 
to a canonically equivalent model at C/ > in a magnetic 
field. With the use of a special particle-hole transforma- 
tion, C;-|, (~l)'Ci| and c[y — > c\^, the Hamiltonian (3) 
with [/ < and parameters (/i, B) is mapped onto a for- 
mally identical Hamiltonian with U > Q and the (in gen- 
eral different) parameters (/i', B'). This is especially useful 
for _B = 0, since the new chemical potential is then given 
by^' = f . The main advantage of the mapping is that (for 
symmetric hopping amplitudes) many properties of the 
repulsive model are well known from the literature [16]. 
Furthermore, with the help of a special particle-hole trans- 
formation, superfluid states in the attractive-C/ model are 
mapped onto antiferromagnetic states with a staggered 
magnetization in the a; — y-plane in the repulsive-C/ model, 
while charge density wave states (CDW) are mapped onto 
antiferromagnetic states with a staggered magnetization 
in z-direction. 



2.2 Diagonalization of the Hamiltonian at 
Hartree-Fock level 

In order to analyze superfluid states at weak coupling, we 
diagonalize the Hamiltonian (2) within the Hartree-Fock 
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approximation by decoupling the interaction term n-^^n-^ 
in the usual manner: 



(4) 



where Ua- = (j^ig.) ^md A = (c-^^c-^^) are assumed to be 
translationally invariant thermal averages. The resulting 
mean-field Hamiltonian can as usual be diagonalized by a 
Bogoliubov-transformation. One finds 



-^fu{\A\' 



'k(T 



(5) 



3.1 Particle numbers at fixed A 

Integral equations as a mapping. In order to determine 
the number of solutions of the coupled equations (7a) 
- (7c), we first consider Equations (7a) and (7b) sepa- 
rately at fixed order parameter A. We show that, at fixed 
A, there is a unique solution of (7a) and (7b) by using 
Banach's fixed point theorem. As vectors we use := 
(riif , nii)"^, where i denotes different points in the para- 
meter space, i.e., different values of {U, t, fi, B}. The inte- 
grations of (7a) and (7b) can be interpreted as a mapping 
X from the parameter space onto itself. Furthermore, in 
order to apply Banach's theorem, we consider the maxi- 
mum metric \\ni — nj||oo = niax£r{|ni,o- — JT-j.o-l}- To prove 



kcr 



E^{e) = sig^{E{e))^me)VlP\W + 

Uini - n|) 



+ aUt^-t^)e + B + 



(6a) 
(6b) 



where a is alternatively interpreted as t,i or ± in (6b). 
The three quantities and A can also be determined 

from their definitions as thermal averages. Suppressing all 
details, we just mention as a result that, at zero temper- 
ature, the following three self-consistency equations have 
to be satisfied in the thermodynamic limit: 



1 



UA 



-OO 
OO 

devdie) 

-OO 

devdie) 



\E{e)\sign{-£^{e)) 
v/£^2(e) + C/2|Z\|2 

|ii;(£)|sign(-g^(£)) 
sign(£'(e)) 



^W{e)+TP\A\^ 
X [e{-£^{e))+e{-£^{e))-l] 



(7a) 
(7b) 

(7c) 



where Vdi^) is the d-dimensional density of states of the 
interaction-free Hubbard model. Since particle-hole sym- 
metry is lost if both B ^ and t| ^ ti (see section 2.1), 
we consider in the following only special cases with either 
B = or i| = ti. 



3 Properties of the self-consistency equations 
for symmetric hopping (t^ = ti = t) 

Since the coupled equations (7a) - (7c) may have more 
than one solution, we show analytically that the number 
of solutions at fixed {[/, i?} is between one and three. 
If more than one solution is found, the one with the lowest 
grand potential is thermodynamically stable. 



uniqueness it is sufficient to show that a number < g < 1 
exists, so that for all n G [0, 1]^: 



where n-^^ is the number operator of the Bogoliubov quasi- 

particles. Moreover, we defined = SiLi cos ki and N is 
the number of lattice sites. The spin-averaged non-superfluid 
particle energies E{e) and the quasiparticle energies £{e) 
are defined as: 



|Z(n + (5n) — In| loo < q | |<5n 



(8) 



holds, where 5n is a small variation in parameter space. 
Hence, it suffices to show that the mapping X is a contrac- 
tion. The solution of (7a) and (7b) is then given by the 
unique fixed point of the equation no = Ing. 



Convergence at small interaction strength. For small 
variations (5n it can be shown that the mapping Z is a 
contraction for all points n in parameter space. In order 
to prove that Equation (8) is indeed satisfied at every 
point n, we use a Taylor-expansion of \T{n„ + 5ncy) — In„\ 
in the parameter 5n and some estimates, whose validity 
depends upon the relative interaction ^ being weak. In 
doing so, we must distinguish two cases, namely that the 
functions ^^(e) and E{e) change sign at the same value 
or at different values of e. In the first case the Taylor- 
expansion yields as a criterion for the validity of (8) with 
< g < 1: 

\U\ 2 

— < ^^-^ (9) 

and in the second case it yields: 



M 
t 



< 



1 



2j^d(en 



(10) 



Here emax is defined as the energy, for which the non- 
interacting density of states is maximal. Since the condi- 
tion (10) is more restrictive than (9), Equation (10) guar- 
antees that the mapping T is a contraction in the whole 
parameter space at sufficiently weak coupling. Note that 
this proof does not work if the density of states diverges at 
a van Hove singularity. For a simple cubic lattice. Equa- 
tion (10) yields ^ < 1.65 for J to be a contraction. 



3.2 Variation of the order parameter A 

In order to solve the coupled Equations (7a) - (7c) we also 
have to discuss the properties of (7c) explicitly. Equation 
(7c) is trivially solved by Z\ = 0, which corresponds to a 
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non-supcrfluid phase. For A ^ Q and symmetric hopping 
ampHtudes it is useful to distinguish between \B\ < \UA\, 
where n-f = rij and Equation (7c) reduces to 



\u\ r°° 1 



U^\A\ 



- 1 



and \B\ > \UA\, where (7c) may be rewritten as: 



M 

2 



v/S2(e) + (72|Z\|2 



(11) 



(12) 



Since, from section 3.1, we know that there is only one 
solution for n| and , we are able to do a one-dimensional 
parameter scan through all possible |Zi [-values in order to 
find the solutions of (7c), using Equations (12) for |Z\| < 
|-^| and (11) for \ A\ > |-^|. Since the integral on the right 
in (11) decreases with increasing \A\ if |Z\| > |-^|, there is 
at most one solution in this interval. Similarly, there is also 
at most one solution for < \A\ < |-^|, since the integral 
on the right in (12) in general increases with |Z\| in this 
case, so that the residual is monotonic on both sides of 
I ^ I . As a consequence, there exist at most three solutions 
for |Z\|, namely the trivial solution |Z\| = and possibly 
two non-trivial solutions on either side of |-^|. Note that 

the modulus of Z\ = ('^ij.'^it) rigorously bounded from 
above by \A\ < i. 



3.3 Possible phases occurring at Hartree-Fock level 

bmce, for |zi| > If |, the magnetization is rigorously zero 
(n-f = ni) and there are at most three non-equivalent 
solutions of the coupled Equations (7a) - (7c), we find 
that (for B ^ 0) there can in principle be up to three 
competing phases: 

— a non-superfluid magnetized phase {A = 0) 

— a superfluid magnetized phase (0 < \A\ < |f |) 

— a superfluid non- magnetized phase (|f | < |Z\|) . 

For i? = it is well known that the superfluid non- 
magnetized solution has the lowest grand potential. Since 
this is not necessarily true for B ^ 0, phase transitions will 
generally be of first order, also at Hartree-Fock level. We 
should also mention that there is no thermodynamically 
stable superfluid magnetized phase at T = and i| ^ ti, 
since that type of solution is always a local maximum of 
the grand potential which can be shown with elementary 
mathematical methods. 



4 Numerical results for symmetric hopping 

Numerical method. With the knowledge of the proper- 
ties of the self-consistency equations (7a) - (7c), estab- 
lished in section 3, it is possible to do accurate numer- 
ical calculations for the phases occurring in the ground 
state. Specifically, we first determine all possible solutions 



of Equations (7a) - (7c), and then identify the thermo- 
dynamically stable solution by comparing the respective 
numerical values of the grand potentials. The various so- 
lutions of the selfconsistency equations are determined by 
performing a scan over |Z\|-values between and ^ in or- 
der to find the roots of (7c). The parameters n-f and nj^ are 
varied self-consistently at each step by successive approx- 
imation; this method is guaranteed to work on account of 
Banach's fixed point theorem (see section 3.1). The phase 
diagram of the homogeneous system is presented in Figure 
1. The parameters n, M and A are presented as a function 
of n and B at fixed U and t. At B = the system is always 
superfluid and the parameters n, M and A do not depend 
on B at flxed fi until a critical magnetic fleld Be, where 
superfluidity breks down and magnetization occurs. Note 
that this first order phase transition generally takes place 
at \Bc\ < |J7Z\_B=o|- The superfluid order parameter is of 
course maximal at half filling {p, 



2 . 



Numerical results within the LDA. In order to formulate 
predictions for experiments on ultracold quantum gases 
we also have to treat the parabolic potential arising from 
the magneto-optical trap. In Hubbard model language this 
potential may be described as: 



^^Trap = ^(i • 12 • i)?!;^ 



(13) 



where 12 is a real, symmetric, positive definite matrix. 
We assume the trapping potential to be spin-independent. 
The trapping potential is treated within the local den- 
sity approximation (LDA), meaning that local quantities 
at site i are determined from the Hartree-Fock selfcon- 
sistency equations corresponding to the effective chemical 
potential /i(i) = fiQ — {i ■ f2 ■ i); the global chemical poten- 
tial fiQ is then tuned so as to reproduce the desired total 
number of particles in the trap. 

The numerical results, obtained from the Hartree-Fock- 
LDA approximation, are presented in Figures 2-4, which 
show the spatial distribution of the occupation number 
n = n^' -\- ni, the spatial distribution of polarization M = 
n| —ni and the spatial distribution of the superfluid order 
parameter's absolute value respectively. In accordance 
with the ellipsoidal geometry of the trap, we flnd a fasci- 
nating and intuitively extremely plausible shell structure 
for the spatial distribution of the various densities and the 
local order parameter. In Figure 2 there is a normal phase 
in the center of the trap surrounded by a superfluid non- 
magnetized shell and a normal phase in the outer region 
of the trap. Apparently, the normal phase in the center 
of the trap is stabilized by the Zeeman-like magnetic fleld 
and reacts very sensitively to changes in this fleld. This can 
be seen from Figure 3, where the (absolute value of the) 
magnetic field is slightly lower than in Figure 2; also the 
global particle number (chemical potential (Iq) is slightly 
different. These changes suffice to cause the disappearance 
of the non-superfluid core. By enhancing the stiffness of 
the trap potential in the x- and z-directions or slightly 
decreasing the interaction strength \U\, it is possible to 
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suppress superfluidity alltogether: There is no superfluid 
phase in Figure 4. 



Experimental consequences. In experiments, the super- 
fluid order parameter A cannot be observed directly. On 
the basis of the results, presented above, it seems that 
little information of use concerning superfluidity can be 
deduced from the total particle number ri| However, 
the other observable considered above, i.e., the local po- 
larization M = n| — rij, might in fact be quite helpful 
in identifying the shell-structure in the trap as an exper- 
imental signature for superfluidity in unbalanced Fermi- 
mixtures. The observable M ~ — can be resolved 
spatially with the use of in-situ imaging [1,2]. 



Comparison with optical lattice-free systems. Systems 
without optical lattices have been intensely discussed in 
the literature, both from a theoretical [3,4] as well as from 
an experimental point of view [1,2]. Important similari- 
ties between systems with and without optical lattices are, 
flrst, the occurrence of superfluidity below a critical tem- 
perature Tc and secondly, the observation of a shell struc- 
ture containing normal and superfluid rings. Moreover, 
there are also several differences between these classes of 
systems, which will now be highlightned. 

In optical lattice-free systems it is found, in contrast 
to our results, that, if superfluidity occurs, the atoms in 
the center of the trap are always superfluid. This differ- 
ence arises from the fact that, in the one-band situation 
discussed in this paper, the quasimomentum is bounded 
from above by ki < 2tt, whereas in optical lattice-free sys- 
tems there is no upper bound for the components of the 
momentum vector. We have shown that superfluid mag- 
netized phases occurring in optical lattice-free systems as 
demonstrated, e.g., by Parish et al. [17,18] do not occur 
in systems, considered in this work (see section 3.3). 

Deformations of the shell structure in the optical lattice- 
free system, discussed by H. Stoof in [3], seem to arise from 
a surface tension between a superfluid core and a normal 
state, which is not discussed in our work, where we con- 
sider systems described by discrete lattice vectors. 

Moreover, in recent publications [19] FFLO phases in 
lattice-free systems are discussed. At present it seems not 
possible to combine this formalism, used for spatially homo- 
geneous (or periodic) potentials, with the LDA-approach, 
used in this paper for spatially inhomogeneous systems, 
since the FFLO-ground state is spatially inhomogeneous 
by itself. For this reason, and since FFLO phases were not 
yet observed experimentally, we disregard FFLO phases in 
this paper. 



Trapping potential-free systems. In recent literature [20], 
trapping potential-free systems have been discussed in the 
weak as well as in the strong coupling regime. In the weak 
coupling regime we obtain the phase diagram in the grand 
canonical ensemble, presented in Figure 1. Since we flnd 



a first order phase transition between a superfluid non- 
magnetized phase and a normal phase, we can formulate 
predictions for systems with flxed particle numbers n-^ and 
ni instead of fixed conjugate variables ^ and B: For non- 
magnetized systems n-f = the system is completely 
superfluid. For a population imbalance, n| ^ nj^, the fol- 
lowing case differentiation has to be done: 

1. If in the normal phase (i.e. assuming A = 0) the par- 
ticle numbers n-^ and can be realized (within the 
grand canonical ensemble) by choosing definite values 
fi and B for the chemical potential and the magnetic 
field, respectively, then this normal phase is also the 
ground state. 

2. If in the normal phase there exists no such combination 
of (/i, i3)-values, then phase separation occurs in the 
ground state. The superfluid and the normal fraction 
have to be determined by a Maxwell construction. 

As mentioned in section 4 a superfluid magnetized phase 
does not exist, at least not in the weak coupling limit. 

5 Results for asymmetric hopping 

In this section we discuss the properties of the Hartree- 
Fock Hamiltonian for superfluid states with asymmetric 
hopping (i| 7^ ^x)- -F'or clarity and convenience, we focus 
on translationally invariant solutions (i? = 0) for _B = 0. 

5.1 Charge density wave states and the repulsive U 
model 

As emphasized before, we are interested in (and, hence, fo- 
cussing on) superfluid phases in the negative-C/ Hubbard 
model, disregarding possible competing phases, which seem 
to be irrelevant for experiments on ultracold quantum 
gases (see for example [11, 12,21-24]). The issue of possible 
competing phases becomes particularly relevant for B = 
and half-filling, since in this case charge density wave 
(CDW) states arc known to be degenerate with s-wavc 
superfluidity if t^ = t^. This is immediately clear from a 
special particle-hole transformation from the negative-?/ 
to the positive-?/ Hubbard model, both at half-flUing and 
-6 = 0, since under this transformation s-wave and CDW 
states are mapped onto the various directions of the stag- 
gered magnetization, which are energetically degenerate 
due to the Hubbard model's S'[/(2)-symmetry in the spin 
sector. Away from half-filling, however, s-wave superfluid- 
ity is stabilized with respect to the CDW phase. 

For moderately asymmetric hopping (i| — ^x) t^*^ sit- 
uation is similar: Exactly at half-filling the CDW phase is 
energetically even somewhat lower than s-wave superflu- 
idity, so that one expects the CDW state to dominate in 
the phase diagram in a narrow but finite strip around half- 
filling. We conclude that the negative-C/ Hubbard model, 
considered here, predicts phases near half-filling, which at 
present seem to be of little experimental interest. For this 
reason we concentrate in the following on band fillings, 
which deviate signiflcantly from half-flUing. so that the 
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superfluid phases, considered here, can safely be assumed 
to be stable [16]. 



5.2 Hartree-Fock density of states 

The absence of symmetry in the hopping amplitudes {t-\ ^ 
ti) clearly also influences the (Hartree-Fock) density of 
states, which is defined as: 

= r deva{e)5{E - E,{e)) . (14) 

The asymmetric hopping now causes additional singulari- 
ties (jumps) in the spin-dependent densities of states. If 
we assume < t| < i|, the following changes occur: 

— The inverse square-root divergence of the density of 
states at the superfluid gap in Vd,'\ is replaced by a 
discontinuity (jump). 

— The DOS Vd^i retains its inverse square-root divergence 
but now displays two additional discontinuities within 
the band (one below and one above the gap). 

— The superfluid gap in Vd,i for i| ^ is smaller than 
the corresponding value 2|[/Z\| for = ^x, while the 
gap remains unchanged in Vd.]- 

The normalization of Vd^a is obviously not affected. 

The densities of states for both spin species, as deter- 
mined from Equations (14) and (6b), are shown in Fig- 
ures 5 and 6 for the parameter values {U = —5, t| = 
3, ti — 0.5, /i = —2, rij = 0.7, ni = 0.3} and a three- 
dimensional simple cubic lattice. These parameter values 
are only used to illustrate the mentioned effects on the 
DOS; they do not necessarily represent a solution of the 
self-consistency equations (7a) - (7c). The spin-dependent 
size of the gap, as seen in Figures 5 and 6, can of course 
also be observed experimentally. 



0.12- 




Fig. 5. Hartree-Fock DOS of the spin species with the greater 
hopping amplitude. Instead of a square-root singularity there 
is a jump at the border of the superfluid gap. 
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Fig. 6. Hartree-Fock DOS of the spin-species with the lower 
hopping-amplitude. Note that an additional jump appears 
within the band, and that the superfluid gap is smaller than 
2\UA\. 

5.3 Phases in the ground state 

In this section we present the results of the numerical so- 
lutions of (7a) - (7c) for translationally invariant simple 
cubic lattices at fixed particle numbers n = -I- and 
B = 0. The proof presented in section 3.1 is also valid for 
asymmetric hopping (i| ^ t|), albeit for smaller values 

of the dimensionless interaction parameter — ^'^L -, than 
for systems with symmetric hopping. This is due to the 
fact that the quasiparticle energy function f (e) in (6b) is 
not a monotonic function in e for the case of asymmet- 
ric hopping. Numerical calculations show that, for i? = 
and t| 7^ ij,, Equation (7c) has at most one solution with 
/I 7^ 0. If two solutions (one with Z\ = and one with 
Z\ 7^ 0) exist, we find that the superfluid phase always 
minimizes the free energy. 

In Figures 7 and 8 we illustrate the behavior of the 
chemical potential u, the magnetization m = — = , 

and the absolute value of the superfluid order parameter 
|Z\| as a function of the interaction strength U at differ- 
ent fixed occupation numbers n = n-f -I- . We find that 
superfluid solutions exist only for attractive interaction 
strengths above a certain minimal value, —U > ?7o > 0. 
In contrast to systems with spin-independent hopping, we 
now also find thermodynamically stable magnetized su- 
perfluid solutions. Note that the trasitions occurring here 
are in general of second order. 



To summarize, we analyzed a generalized attractive-[/ Hub- 
bard model with (possibly spin-dependent) hopping am- 
plitudes and a Zeeman-like magnetic field, which is rele- 
vant for ultracold quantum gases. In view of the spatial 
inhomogeneity due to the presence of a trap, we stud- 
ied the model at Hartree-Fock level. For the special case 
of symmetric hopping we analytically demonstrated the 
uniqueness of the solution of Equations (7a) and (7b) at 
fixed order parameter A. Hence we are certain to have 
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detected all solutions of (7a) - (7c). We found that phase 
transitions at i? ^ are generally of first order. Within 
the local density approximation we showed that phase 
separation with an interesting shell structure occur if a 
parabolic potential is added to the Hamiltonian (2). For 
general spin-dependent hopping and i? = we calculated 
the (spin-dependent) density of states. Interestingly, we 
found that, sufficiently far away from half filling, a super- 
fluid magnetized phase occurs, which was absent in the 
phase diagram for spin-independent hopping. 

Finally, we discuss possible extensions of our results. In 
this paper, we restricted consideration to weak-coupling 
solutions at Hartree-Fock-LDA level. For our purposes, 
this restriction is meaningful, since the Hartree-Fock-LDA 
approximation is expected to yield qualitatively correct re- 
sults in the weak coupling regime, provided spatial gradi- 
ents are not too large; moreover, any other method would 
be virtually impracticable in a spatially inhomogcneous 
system with a magneto-optical trap. Nevertheless, with an 
eye on the future, it seems worthwhile to review possible 
alternatives. 

First of all, it would probably be possible to drop the 
LDA-approximation and calculate the spatially inhomo- 
gcneous Hartree-Fock solution for an optical lattice in a 
parabolical trap. However, even for a two-dimensional lat- 
tice without spin-asymmetries (5 = and t| = i|), this is 
known to be computationally fairly expensive for realistic 
lattice sizes (with typically 10'^ particles) [25]. 

Secondly, it would of course be important to go beyond 
the Hartree-Fock approximation and include correlation 
effects. This could be done in principle by taking into ac- 
count the dominant quantum fluctuations, which arc de- 
scribed by second order Feynman diagrams in a pertur- 
bation expansion. In principle, such an expansion could 
be formulated also for spatially inhomogcneous systems 
with symmetry breaking, but the computational expense 
of solving these equations self-consistently for a three- 
dimensional system of reasonable size would be tremen- 
dous. Here we just recall that it is known for translation- 
ally invariant systems [26-28] in three or more dimensions, 
that the leading effect of quantum fluctuations is to renor- 
malize Hartree-Fock results and reduce all energy scales 
(critical temperatures, gaps) by factors of typically three 
to four. For the calculations of the present paper this sug- 
gests that inclusion of second order diagrams would quan- 
titatively but not qualitatively change the Hartree-Fock 
results. 

Clearly it would, for the description of non-perturbative 
and strong-coupling effects, also be desirable to apply sim- 
ulation methods (in particular Quantum Monte Carlo tech- 
niques) to the Hubbard models studied in this paper. 
However, a simulation of a three-dimensional system of 
physically interesting size at sufficiently low temperatures 
seems at present unfeasible. It also seems unlikely that ap- 
proximations like Dynamical Mean Field Theory (DMFT) 
or its cluster extensions, which arc known to be extremely 
helpful in translationally invariant systems in d = 2 and 
c? = 3, are applicable to spatially inhomogcneous systems 
in a trap at ultralow temperatures, where superfluidity 



occurs. However, away from the symmetry-broken low- 
temperature phase, the DMFT has been demonstrated to 
be well applicable, e.g., to the description of spatially in- 
homogeneous Mott metal- insulator transitions [29] . 

For the time being, we conclude that we discovered sev- 
eral fascinating new phenomena in Hubbard-type models 
for asymmetric fcrmionic ultracold quantum gases. Specif- 
ically we found phase separation with an interesting shell 
structure for systems in a trap and also (in a translation- 
ally invariant model) a new superfluid magnetized phase. 
It would be important if the calculations presented here 
could be extended beyond the weak-coupling regime, using 
different methods, along the lines pointed out above. We 
hope and expect that our results can soon be compared 
to experiment. 
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Fig. 1. Pliase diagram of tlie iiomogeneous system a,t T — 0, t — 1 and U — —3. Tiie parameters n, M and A axe sliown as a 
function on and B. The jump in the superfluid order parameter A and the magnetization M indicates the first order phase 
transition. The relative jump of the particle number n is smaller and is therefore not adequate as a signature of superfluidity. 




Fig. 2. Spatial distribution of the parameters n, M = ni — ni and |Zi| in the x-z- plane of a simple cubic lattice at trap 
parameters: {U = -3, t = 1, = 2, B = -0.1, H = Diag(0.05, 0.3, 0.1)}. 




Fig. 3. Spatial distribution of the parameters n, — 7i| — and \A\ in the x-z-plane of a simple cubic lattice at trap 
parameters: {U = -3, f = 1, = 0, B = -0.05, f2 = Diag(0.05, 0.3, 0.1)}. 
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Fig. 4. Spatial distribution of the parameters n and M = ji f — in the x-z-plane of a simple cubic lattice at trap parameters: 
{U — —2, t = 1, ^0 — 1, B = —0.01, n — Diag(0.2, 0.3, 0.1)}. The order parameter \ A\ vanishes in the whole system and is not 
shown. 



^ m A 




Fig. 7. The parameters fi, m and A as a, function of U near half filling n « 0.79 

fJ- m A 




Fig. 8. The parameters fi, m and Zi as a function of U away half filling n ~ 0.18 



